TERMINAL AREA GUIDANCE ALONG CURVED PATHS - A STOCHASTIC CONTROL APPROACH* 


J. E. Quaranta and R. H. Foulkes, Jr. 
Youngstown State University 


ABSTRACT 

In this paper, stochastic control theory is applied to the problem of de- 
signing a digital flight compensator for terminal guidance along a helical 
flight path as a prelude to landing. The development of aircraft, wind, and 
measurement models is discussed along with a control scheme consisting of feed- 
back gains multiplying estimate of the aircraft and wind states obtained from 
a Kalman one-step predictor. Preliminary results are presented which indicate 
that the compensator performs satisfactorily in the presence of both steady 
winds and gusts. 


1. INTRODUCTION 

During the past few decades, the number of problems associated with air- 
port traffic has risen dramatically. Among the more pressing areas of concern 
are high noise levels near airports, fuel conservation, and weather- induced 
delays, diversions, or closures. In order to alleviate some of these problems, 
NASA and the FAA have jointly initiated' a long range research effort, the Term- 
inal Configured Vehicle (TCV) program (ref. 1). Among the objectives of the 
TCV program are increased capability for zero-visibility operation, reduced air 
delays and route time, avoidance of sensitive areas, and reduced noise source 
intensity. These objectives can be met, at least partially, through the devel- 
opment of precise automatic control along steep, curved approach paths. 

A prerequisite to such precise automatic control is the development of im- 
proved ground-based navigation and guidance systems along with improved airborne 
control systems. The ground based improvements in terminal area navigation and 
guidance will be provided by the Microwave Landing System (MLS) . The MLS will 
periodically provide accurate range, elevation, and azimuth information to the 
on-board control ' system. 

The purpose of this paper is to present an application of stochastic con- 
trol theory to the problem of designing an airborne control system that uses 
the MLS data for terminal area guidance of a Boeing 737 along a helical flight 
path as a prelude to landing. First, a system model is presented consisting of 
an aircraft model, wind model, and measurement model. Next, the digital compen- 
sator design is presented. Finally, a digital simulation showing the system 
response using the above compensator is presented. 


2. THE SYSTEM MODELS 

In this chapter the aircraft, wind, and measurement models are developed. 
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2.1 The Aircraft Model 


The notation in this section follows closely that in Etkin C^ef. Thus, 

Fj and Fg denote inertial and Earth reference frames, respectively. If V is the 
mass center velocity of an aircraft, and denote the vector quantity V meas- 
ured with respect to Fj and Fg, respectively. Also is the vector quantity 

expressed in the vertical reference frame (Fy) coordinates. 

The basic nonlinear aircraft model consists of the vector force and moment 
equations, a set of kinematic constraints, inertial^ velocity equa_tions, and a 
set of actuator equations. The force equation is f = ma^, where f is the exter- 
nal force and a^ is the mass center inertial acceleration. Assuming the Earth 
is an inertial system and is locally flat, a^ = (the dot notation means time 
differentiation), so that 'aj.g = + wBg x V^b, where = [p q r]' is the 

angular velocity of Fg and V^g is the mass center inertial velocity, both ex- 
pressed along the body axes. The prime means transpose. 

Let W denote the wind velocity. Also, let F^ be an atmosphere-fixed refer- 
ence frame and V the mass center velocity with respect to Fa- Then^^B “ ^Bw ^w 
+ V?B, where Lg^ transforms wind-axes components into body-axis components 
(ref. 2), Hence, 

^CB "" ''^Bw Vw + Wg) + o)Bg X- (Lg^ + Wg) (2.1-1) 


The external force is f g ^ Ag + Tg + mgg, where Ag is the aerodynamic 
force, Tg is the thrust, and mgg is the gravity force. Hence, the vector force 
equation is 

Ht asw V,, + Wg) + m^g X (Lg^ + Wg) = ^ Ag + ^ Tg + gg (2.1-2) 

The vector moment equation is G = h, where G is the external moment and h 
is _the_sngular momentum. Assuming 1=0 and neglected and elastic components 
of h, hg = I cJ®g, where I denotes the body-axes moments of inertia. Thus, 

Gg = I mBg + w®g X I m^B (2.1-3) 

The external moment Gg includes the effects of the gust angular velocities 
Pg, qg, and rg. These effects show up explicitly in the perturbation model. 

The scalar components of (2.1-2) and (2.1-3) involve both the body-axes 
Euler angles (f>, 6, and ip (bank, pitch, and heading angles, respectively) and 
the body-axes rates p, q, and r. The two sets of variables are related by 

<j) = p + q sin (J) tan 0 + r cos tan 0 

0 = q cos (j) - r sin <j> (2.1-4) 

Tp = (q sin (j) + r cos cji) sec 0 


The Earth position is described in cylindrical coordinates with origin at 
the helix center at ground level. The position of the origin is assumed to be 
known with respect to the MLS origin. The rates of the Earth position coordin- 
ates (helix radius R, helix angle v, and altitude h) are given by 


R = V cos Y cos (ijj - v) + Wx cos v + Wy sin v 
V Wx Wy 

V = g cos Y sin (ijj^ - v) - -g— sin v + ^ cos v (2.1-5) 

h = V sin Y “ W 2 
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where y and are wind-axes elevation and heading angles and Wx, Wy, and Wz 
are coordinates of W in the reference frame 

The remaining equations result from modeling the thrust throttle and sta- 
bilizer actuator systems. The thrust- throttle relation is modeled as a first- 
order lag with a time constant of 0.5 second. Because this relation is a lin- 
earization about nominal values, the equation is given with the development of 
the perturbation model. In addition, throttle and stabilizer rates are command- 
ed inputs, so that throttle and stabilizer positions are state variables. 

The nonlinear aircraft model consists of the force equation (2.1-2), the 
moment equation (2.1-3), the kinematic constraints (2.1-4), the inertial velo- 
city equations (2.1-5), and the actuator relations. The model can be written 
in usual state variable form as a single, nonlinear vector equation 

X = f(X, U, W, W) (2.1-6) 

where X_= [VBapqr<j)0ij)RvhT7r4]' is the total state vector, 

U = [tt -5 6e dr 6a ^sp]' ig the commanded input vector, W = [ug Vg Wg pg qg rg 
Wr Wt] ' is the wind vector, and where 6 is the sideslip angle, a is the angle 
of attack, T is the thrust, tt is throttle setting, 4 is stabilizer, 6g is the 
elevator, gj- is the rudder, ,5a is the aileron, ^gp is the spoiler, Ug, Vg, and 
Wg are translational gust velocities, pg, qg, and rg are rotational gust velo- 
cities, and and W'p are radial and tangential components of the steady wind, 
respectively. 

The perturbation model consists of the first-order terms in a Taylor 
series expansion of (2.1-6) about a descending helical equilibrium. The equi- 
librium was determined under a zero wind condition with the aircraft flying a 
"truly banked" turn (ref. 2) for an airspeed of 64 m/sec (120 knots), bank angle 
of 15° and angle of elevation of using data for the Boeing 737 from the TCV 

program. The coefficients in the perturbation model were computed by evaluating 
the appropriate partial derivatives at the equilibrium. The coefficients of 
pg, qg, and rg in the moment equations were computed by evaluating the partials 
of the aerodynamic terms in (2.1-3) with respect to p, q, and r, respectively. 
Finally, the thrust- throttle relation is 6T = - 0.5 6T + 313.33 6ir, where the 
thrust is in pounds, the throttle setting in degrees, and the "6" indicates a 
perturbation value. The perturbation model in usual state variable fom is 

X = Ax + Bu + Dqiw + Dq2 w (2.1-7) 

where x, u, and w are the perturbation counterparts, of the total vectors X, U, 
and W in (2 . 1-6) . 

2 . 2 The Wind Model 


As seen in section 2.1, the wind vector in the aircraft model consists of 
three translational gust velocities ug, Vg, and Wg; three rotational gust velo- 
cities Pg, qg, and rg; and two steady wind components Wr and Wx. The gust velo- 
cities, all of which are components along the body-axes, are modeled as having 
the Dryden spectra and are produced for simulation and filter design purposes 
by a linear system processing white noise. As an example of the linear system 
design, consider the gust velocity Ug, normalized by the equilibrium airspeed 
Ve. The power density spectrum of the normalized Ug is “ (^L^Ou^/Vq ) 

/ (1 + (LuOj/Ve^) , where cTu is the rms gust velocity, is a turbulence scale 
factor, and m is the frequency variable in rad/sec. Now, if a linear system 
with transfer function H(joj) = 1/(1 + jo) Lu/Ve) is subjected to a white noise 
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input with variance = ( 2 Lu 0 u 2 /Ve 3 ) , the output is a random process with the 
spectrum (f)u (w) (ref. 3). A system with the required transfer is described in 
state-variable form by x^l = (-Vg/Ly) + (Vg/Lu) = x^]^, where is a 

mean zero white noise process with variance ( 2 Luau^/Vg^) , x^]^ is a state vari- 
able, and Wq^ is the output having the required spectrum. The remaining gust 
velocities are generated in a similar manner. 

As indicated in the previous section, and are the radial and tangen- 
tial components of the steady wind, which are related to the north wind and 
east wind Wg by the spiral angle v: 


Wr = Wjvj cos V + Wg sin v ^2 2 - 1 ) 

Wx = W{(j sin V - Wg cos v 

Thus, Wr = - V Wf and W-p = v Wr. For the simulation purposes, north and east 
winds are selected and Wr and W'p are computed using equations (2.2-1). For the 
filter design, because a cpnstant ^ coefficient wind model is desired and the 
equilibrium wind is zero, Wr and W-p are approximated by 


Wr = - Ve Wt 
Wt = Vg Wr 


( 2 . 2 - 2 ) 


where Vg is the equilibrium spiral angle rate, which is constant and equal to 
Vg cos Ye/^e- 

Putting together the gust system equations and the steady wind equations 
( 2 . 2 - 2 ) yields a time- Invariant, linear wind model of the form x^ = A^ x^ + Bw^ 
and w = x^, where ? is the white noise vector generating the gusts, x^ is the 

state vector of the wind model, and w is the wind vector used in the perturba- 
tion^ model. In order to use x^^ in the perturbation model equations, w = c^ x^^ 
and w = Av 7 Xw + Cw Bw ^ are substituted into equation (2.1-7) to give a com- 
bined aircraf t/wlnd model: 


X = Ax + Bu + Dq x^^ + 5 

Xw = Aw ’^w + ^ ^ 

where Dq = C„ + Dq 2 \ and Dj = Dq2 B^ 


2.3 The Measurement Model 


Measurements available for control purposes consist of the MLS data (range, 
azimuth, and elevation) and a number of on-board sensor readings. The total 
measurement vector is Y = [^ Az p q r (f* © V h|j hb xg yg zg] ' , where -Y, Az, 
EZ are the MLS data; p, q, r are angular velocities from rate gyros; <!>, 9, ip are 
bank angle, pitch, and heading from position gyros; V is an airspeed indicator 
reading, h^ and h|j are barometric altimeter and vertical speed indicator read- 
ings; and xg, yg, zg are body-mounted accelerometer readings. 

The total measurements are computed for simulation purposes by computing 
the total state variables as equilibrium values plus increments and expressing 
the measurements in terms of the states. In order to compute the MLS data, it 
is assumed that the A, Az, and EZ are measured with respect to a common origin 
and that the helix center is known with respect to that origin. If the ground 
coordinates of the helix center with respect the MLS origin are (Xq, yo ) 5 then 


1770 



fl = 


° 


-1 


Az = tan 


EZ. - tan ^ 


+ R cos v)^ 4- (yq + R sin v)2 + 
Yq + R sin 


[ 


o 


V 

+ R cos V 


r 

L (Xc 


) 

) 


■1 


+ R cos v)2 + (yq + R sin v)^ 
where R, v, and h are coordinates 10, 11, and 12 of the total state vector. 

The fourth through eleventh measurements are the same as total states. 

Also h]j = Vg sin Ye ^12’ '^^ere the derivative xj2 computed from equation 
2.2-3). FinallY, expressions for the accelerometer readings are obtained bY 
writing out the scalar components of the acceleration from equation (2.1-1). 

In the simulation, the measurements are generated using the above relations 
along with random noise effects. Except for the airspeed and vertical speed 
indicators, the noise is an additive, white, mean zero Gaussian process with 
standard deviation as shown in Table 1 (ref. 4 and 5). The airspeed and verti- 
cal speed indicator noises are multiplicative, where the indicated measurement 
is obtained bY multiplYing the actual measurement bY a. normal random variable 
of mean 1 and standard deviation as given in Table 1. 

The incremental measurements to be used with the perturbation model are 
calculated bY subtracting the equilibrium measurement values from the total 
measurements discussed above except for the first three incremental measurements. 
For these first three, total helix radius R, helix angle v, and altitude h are 
computed from range, azimuth, and elevation using the following equations: 


R = 


= (A. cos E£ cos Az - Xg)^ + (A. cos EZ. sin Az 

fi cos EZ, sin Az - Y 


Yo)2 


= tan 


-1 


Y 

h = A. sin EZ. 




fi cos EZ sin Az - 


X 


0 


Then the equilibrium radius, angle, and altitude are subtracted to generate the 
incremental measurements. 

The fifteen incremental measurements are linear functions of x and x^^. 
Hence, the incremental measurement vector y can be written as y = Cx + Cw ^w+ v, 
where v is a noise term whose coordinates are assumed to be white, mean zero 
Gaussian processes with standard deviations that are the same as for the total 
measurements except for the noise terms in yi> Y2’ Y3’ YlO’ Yl2 (see Table 

1 ). 


3. THE CONTROL SYSTEM DESIGN 

Using the aircraft, wind, and measurement models presented in the previous 
chapter, the total system model is 

x(t) = Ax(t) + Bu(t) + Dgw(t) + Di^(t) (3-1) 

w(t) = A„w(t) + B„4(t) (3-2) 

y(t) + Cx(t) + C„w(t) + v(t) (3-3) 

where x, u, w and y are the state, control, wind disturbance, and measurement 
vectors resp., v(t) is a white, Gaussian vector of measurement noise and 4(t) is 
a white, Gaussian noise vector that drives the wind system and corrupts the air- 
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craft system. The matrices A, B, D^, A^, B^, C and are time- invariant 

with appropriate dimensions. 

The problem of designing a feedback compensator is posed as the usual lin- 
ear, stochastic regulator problem. For the regulator problem, a quadratic cost 
functional of the form 

J = -7 E {/ ^ [x'(t) Qx(t) + u'(t) Ru(t)]dt} (3-4) 

to 

is used, where E is the expectation operator, the prime denotes the transpose, 
and both Q and R are positive definite, time- invariant weighting matrices. The 
problem can now be stated as follows. Given the linear system of equations 
(3-1 to 3-3), find a control u such that the cost functional J of equation (3-4) 
is minimized. 

The first step in solving this problem is to transform the system of equa- 
tions (3-1 to 3-4) into their discrete-time equivalents. This is done for 
several reasons. First, a digital compensator is desired since the on-board 
computer is digital and any control algorithm must be compatible with a digital 
system. Secondly, the MLS data is only provided periodically. Therefore, the 
measurement system is inherently a discrete-time one. Finally, a digital simu- 
lation is used to test each design. Therefore, the discrete-time equivalent 
difference equations make the simulation very easy to implement on a digital 
computer. 

The equivalents are obtained by integrating the system differential equa- 
tions and cost functional over each sampling period (ref. 6 ). 

If we restrict u(t) to be constant over the sampling period, the resultant 
discrete-time equations are 


^k+1 = <f>Xk + ^2 u^ + Ck (3-5) 

Wk+1 = <t>w + Tij^ (3-6) 

where the subscript k denotes the kth sample and | ' is a mean zero, 

white noise sequence whose variance depends on the variance of the continuous 
noise vector ^(t). 

The discrete-time equivalent measurement equation can be obtained directly 
from the continuous time equation: 


y = Cxk + C„wk + vk (3-7) 

Finally, the discrete-time equivalent of the cost functional (3-4) can be 
written as a sum of n integrals. The resultant expression for the cost func- 
tional becomes : 


" t ^ ^k=o ^k+1 Q%+1 + 2x^+1 Nwj ^+1 + 2x^ Muk + u^ Ruk> (3-8) 


Note that, since the original system and cost matrices are time- invariant , 
so are the discrete-time system and cost matrices. 

The problem can now be restated as follows. Find a control sequence uk 
which minimizes the cost functional J in equation (3-8) subject to the con- 
straints that the state equations (3-5, 6 , 7) must be satisfied and that (uk) 
must explicitly depend only on the past measurements yj^ = {yg, yj^, . . . , yk-l^> 
where vk is a zero-mean, Gaussian, white noise sequence independent of 

The above system of equations (3-5, 6 , 7, 8 ) can be augmented to obtain a 
form very similar to the discrete linear quadratic Gaussian stochastic control 
problem. However, if the normal method of solution is applied, an important 
difficulty surfaces. The total system may be unstable and uncontrollable due 
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to an unstable wind system. Therefore, if the augmented system is solved with 
an unstable wind system, the solution to the Riccati equation diverges due to 
the presence of unstable and uncontrollable poles. But, under certain con- 
ditions the gains will be bounded. 

It can be shown (ref. 6 ) that the solution to the stochastic optimal con- 
trol problem described previously exists and is given by; 

= - Hk^k - H^k^k 

% = Hwk = ^wk> \ = R + ri 


Gk = Pk‘*> + M’ , = r- + Pj^r2) 

Pk-1 = <i>'Pk'<> + Q - Gi Gi,, Pn = Q 


Pwk-1 = (<j>-riR-ir^) ’ (p^^<f.„ + Pkr2) + n, = n 

where Xj^ and Wj^ are one-step predicted estimates of Xj^ and wj^ given by; 
% = E {x^|yj^_3^} and \ = E 


It should be noted that the above gain equations remain valid for any i-step 
predicted estimate (where i = 0 represents a filtered estimate) . The one-step 
prediction- was used here in order to account for computational delays present 
in the on-board computer. 

Also, the equations above are of a recursive nature. Therefore, at each 
sampling instant a new optimal gain is calculated. To implement this would 
require storing all the intermediate values of each gain matrix. This, in turn 
would require a greater amount of storage than is normally available for small, 
on-board computers. For these reasons, a suboptimal design was used consisting 
of only the steady state gains obtained when the index on the recursive rela- 
tions tends to infinity. 

It should be emphasized that the optimal control for a system with distur- 
bances consists of two parts. The first part feeds back state estimates multi- 
plied by an optimal gain H^. This gain is exactly the same as would have been 
calculated with no disturbance present. The second term feeds back the distur- 
bance estimates multiplied by a gain which depends on the disturbance. 

The next step is to obtain the state and wind estimate. This is accom- 
plished by first augmenting the discrete-time equation^ (3-5, 6, 7): 


r^k+ri r 

ft] • [f 

0 • ft] 

r 

Yk [G}C^j] 

0 



+ Xk 


X]^ is a white noise cerm representing modeling error. 

Given the past measurements y^^, it can be shown (ref. 7 ) that the one-step 
predicted estimates of Xk+i and Wk+q are given by the equations below; 

2 ' 


"^k+ll ^ 
-^k+l_ 


where = 




= E {x,}. 
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^k 
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^-"k. 

0 , and 
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where 


= 



tcSc^]’ 


([cic^] [c!c^]' 


0k) 


-1 



r/v- -| 

^k 

f 

} 

?k 


^k 


^k 

LwkJ 



-'^k- 


-'^k- 


A 

LwkJ 


is the estimation 


error, and Zj^ can be found by solving the Riccati type equation: 


where 


and 



([clc^] Sk [c‘c„]' 


= E { 


} = E { 



"5k ' 

-0k- 

-0k- 


Isl 


1 

1 

o 

LikJ 


L%J 


+ Qr) ^[c! c„] Zi,} 


} + E {x^ x{} 


} = 0 for all k 


j =0, 1, 2, 

It must be noted that, as with the control equations, these equations are 
recursive. For the same reasons discussed previously, a suboptimal predictor 
was implemented using only the steady state solutions to the above equations. 

Now that the optimal controls have been defined for specific cost matrices 
the main concern becomes testing to see if the chosen cost matrices lead to an 
acceptable system response or if they must be modified to achieve this goal. 

The digital simulation described in the next section will provide the final 
step of the design procedure satisfying the above testing and modification 
requirements . 


4. SIMULATION RESULTS 

The main objective considered in generating these results was to design 
the control system so that the aircraft position was kept close to the equilib- 
rium in the presence of various steady wind magnitudes. Thus, the objective was 
to keep the helix radius, altitude, and airspeed perturbations small while 
allowing some of the other perturbations, such as the altitude variables, to 
become relatively large. 

The first set of results (Table 2) shows how the control system reacts to 
increasing steady winds. Three simulation runs were made with steady wind velo- 
cities of 3.05 m/sec (10 ft/sec), 6.1 m/sec (20 ft/sec), and 12.2 m/sec (40 
ft/sec). The rms gust velocity in all runs was 0.61 m/sec (2 ft/sec). The 
results in Table 2 indicate that, as the steady wind increases, the maximum 
deviations in the bank angle, heading angle, and helix angle also increase pro- 
portionally, while the maximum deviations in helix radius, altitude, and air- 
speed increase at a much smaller rate. 

The results in Table 3 show the importance of the control gain matrix 
in controlling the aircraft in the presence of winds. Two runs were made with 
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12.2 m/sec (40 ft/sec) winds, one with used in the control scheme and one 
without H„. As the results indicate, the term in the control calculations 
has a large effect in providing satisfactory control in the presence of winds. 


5. CONCLUDING REMARKS 

A linear, time- invariant perturbation model of an aircraft flying a des- 
cending helix in the presence of winds was developed. An automatic control 
system was designed through an application of the linear-quadratic-Gaussian 
discrete- time regulator theory. 

Using a performance criteria including maintaining a circular ground 
track of correct radius as well as maintaining proper altitude and airspeed, 
preliminary simulation results have shown that the control system can perform 
satisfactorily in the presence of steady winds without excessively large devi- 
ations in other variables, such as bank and pitch angles. 
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TABLE 1. MEASUREMENT NOISE STANDARD DEVIATIONS 


MEASUREMENT 

NOISE STANDARD DEVIATION 



Range 

.305 m (1 foot) 

Azimuth 

0.41 X 10“^ rad 

Elevation 

0.61 X 10 “^ rad 

Helix Radius 

1.52 m (5 feet) 

Helix Angle 

10~3 rad 

Altitude 

3.05 m (10 feet) 

Rate gyros (p, q, r) 

0.1 deg/sec 

Bank gyro 

0.5 deg 

Pitch gyro 

0.15 deg 

Heading gyro 

1 . deg 

Airspeed (total) 

.02 

Airspeed (incremental) 

.02 V 3 

Barometric altimeter 

8.38 m (27.5 feet) 

Vertical speed Indicator (total) 

.05 

Vertical speed indicator 
(incremental) 

.05 hg^ 

Accelerometer (xg, Zg) 

4.91 cm/sec^ (.161 ft/sec^) 

Accelerometer (yg) 

.491 cm/sec^ (.016 ft/sec^) 
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TABLE 2. EFFECT OF STEADY WINDS 


Steady Wind Magnitude, m/sec (ft/sec) 

3.05 (10) 

6.1 (20) 

12.2 (40) 

Maximum Perturbation Magnitude 




Bank angle, deg 

2.0 

3. 1 

6.7 

Heading angle, deg 

9.1 

16.3 

31.4 

Helix angle, deg 

5.2 

9.7 

19.5 

Airspeed, m/sec (ft/sec) 

1.7 (5.6) 

1.4 (4.5) 

2.0 (6.6) 

Helix radius, m (ft) 

5.6 (18.3) 

7.3 (24.0) 

11.0 (36.0) 

Altitude, m (ft) 

7.3 (24.0) 

8.1 (26.7) 

10.0 (32.8) 


TABLE 3. EFFECT OF CONTROL GAINS 


Maximum Perturbation Magnitude 

with H^ 

without H^ 

Airspeed, m/sec (ft/sec) 

2.1 

(6.8) 

3.2 

(10.4) 

Helix radius, m (ft) 

11.3 

(37.0) 

81.7 

(268.1) 

Altitude, m (ft) 

10.3 

(33.8) 

16.0 

(52.4) 



















